% KurmannSims_Part1.m
%
% Code to replicate Tables 1 and 2 of main text as well as Table 1 and Figures 1-3 of Appendix of
%
%   Kurmann, A. and E. Sims (2019). "Revisions in Utilization-Adjusted TFP and Robust Identification of News Shocks." Review of Economics and Statistics (forthcoming).
%
% Written by Andre Kurmann (Drexel University) and Eric Sims (Notre Dame and NBER) except for
%   - NBERbc.m (plots time series with NBER recession bars)
%   - bwfilter.m (Stock and Watson's bi-weight filter as used by Fernald)
%   - cf_bpass_default.m (Christiano-Fitzgerald bandpass filter as used by Fernald)
%
% Last modified: Andre Kurmann, January 10, 2020
%
% Use of code for research purposes is permitted as long as proper reference to source is given. 
%
%------------------------------------------------------------------------------------------------------------------------------------------------

clear all
close all
clc

%% Fernald's different data vintages
tfp_u_data = xlsread('KurmannSims_Part1.xlsx','tfp_util_compare');
u_data = xlsread('KurmannSims_Part1.xlsx','util_compare');
tfp_data = xlsread('KurmannSims_Part1.xlsx','unadjustedtfp_compare');
[T,q] = size(tfp_u_data);
st = 1; %1 for 1947q3, 51 for 1960q1, 147 for 1984q1
et = T; %T for 2007q3, 130 for 1979q4, 
time = [1947.75:0.25:2007.75]; %adjust accordingly

d_tfpu07 = tfp_u_data(st:et,2)*400; %400*log quarter to quarter change  
d_tfpu11 = tfp_u_data(st:et,3)*400;
d_tfpu13 = tfp_u_data(st:et,4)*400;
d_tfpu14 = tfp_u_data(st:et,5)*400;
d_tfpu15 = tfp_u_data(st:et,6)*400;
d_tfpu16 = tfp_u_data(st:et,7)*400;
tfpu07 = tfp_u_data(st:et,8);       %log tfpu (no re-sizing)
tfpu11 = tfp_u_data(st:et,9);
tfpu13 = tfp_u_data(st:et,10);
tfpu14 = tfp_u_data(st:et,11);
tfpu15 = tfp_u_data(st:et,12);
tfpu16 = tfp_u_data(st:et,13);

d_u07 = u_data(st:et,2)*400;    %400*log quarter to quarter change 
d_u11 = u_data(st:et,3)*400;
d_u13 = u_data(st:et,4)*400;
d_u14 = u_data(st:et,5)*400;
d_u15 = u_data(st:et,6)*400;
d_u16 = u_data(st:et,7)*400;
u07 = u_data(st:et,8);          %log u (no re-sizing)
u11 = u_data(st:et,9);
u13 = u_data(st:et,10);
u14 = u_data(st:et,11);
u15 = u_data(st:et,12);
u16 = u_data(st:et,13);

d_tfp07 = tfp_data(:,2)*400;
d_tfp11 = tfp_data(:,3)*400;
d_tfp13 = tfp_data(:,4)*400;
d_tfp14 = tfp_data(:,5)*400;
d_tfp15 = tfp_data(:,6)*400;
d_tfp16 = tfp_data(:,7)*400;
tfp07 = tfp_data(:,8);
tfp11 = tfp_data(:,9);
tfp13 = tfp_data(:,10);
tfp14 = tfp_data(:,11);
tfp15 = tfp_data(:,12);
tfp16 = tfp_data(:,13);

%% Macro data
nipa_data = xlsread('KurmannSims_Part1.xlsx','nipa_clean');
d_y = nipa_data(st:et,8)*400;   %400*log quarter to quarter change
d_c = nipa_data(st:et,9)*400;
d_I = nipa_data(st:et,10)*400;
d_h = nipa_data(st:et,11)*400;
d_h_n = nipa_data(st:et,16)*400;     
y = nipa_data(st:et,2);
c = nipa_data(st:et,3);
I = nipa_data(st:et,4);
h = nipa_data(st:et,5);
h_n = nipa_data(st:et,15);          %log h/n (no re-sizing)
infl = nipa_data(st:et,14)*400;


%% Figure 1 and unconditional moments for Table 1

% Figure 1
figure('Name','Figure 1 of Appendix')
subplot(311);   fh = NBERbc(time, [tfpu07 tfpu16], {'-','-'}, [1.5,1.5], {'k','b'}); 
                title('Utilization-adjusted TFP');
                legend([fh(1),fh(2)],'Adjusted TFP 2007 vintage','Adjusted TFP 2016 vintage')
subplot(312);   fh = NBERbc(time, [tfp07 tfp16], {'-','-'}, [1.5,1.5], {'k','b'}); 
                title('Unadjusted TFP')
                legend([fh(1),fh(2)],'Unadjusted TFP 2007 vintage','Unadjusted TFP 2016 vintage')
subplot(313);   fh = NBERbc(time, [u07-mean(u07)*ones(T,1) u16-mean(u16)*ones(T,1)], {'-','-'}, [1.5,1.5], {'k','b'}); 
                title('Utilization'); 
                legend([fh(1),fh(2)],'Utilization 2007 vintage','Utilization 2016 vintage')

% unconditional moments for Table 1 
format bank
disp('Unconditional moments for Table 1 of main text')

d_tfpu_07 = [mean(d_tfpu07); std(d_tfpu07); corr(d_tfpu07,d_tfpu07)];
d_tfpu_11 = [mean(d_tfpu11); std(d_tfpu11); corr(d_tfpu11,d_tfpu07)];
d_tfpu_13 = [mean(d_tfpu13); std(d_tfpu13); corr(d_tfpu13,d_tfpu07)];
d_tfpu_14 = [mean(d_tfpu14); std(d_tfpu14); corr(d_tfpu14,d_tfpu07)];
d_tfpu_15 = [mean(d_tfpu15); std(d_tfpu15); corr(d_tfpu15,d_tfpu07)];
d_tfpu_16 = [mean(d_tfpu16); std(d_tfpu16); corr(d_tfpu16,d_tfpu07)];
tfpu_momentnames = {'mean'; 'stdev'; 'corr(.,d_tfpu_07)'};
Table_dtfpu_stats = table(d_tfpu_07,d_tfpu_11,d_tfpu_13,d_tfpu_14,d_tfpu_15,d_tfpu_16,'Rownames',tfpu_momentnames)

d_tfp_07 = [mean(d_tfp07); std(d_tfp07); corr(d_tfp07,d_tfp07)];
d_tfp_11 = [mean(d_tfp11); std(d_tfp11); corr(d_tfp11,d_tfp07)];
d_tfp_13 = [mean(d_tfp13); std(d_tfp13); corr(d_tfp13,d_tfp07)];
d_tfp_14 = [mean(d_tfp14); std(d_tfp14); corr(d_tfp14,d_tfp07)];
d_tfp_15 = [mean(d_tfp14); std(d_tfp14); corr(d_tfp15,d_tfp07)];
d_tfp_16 = [mean(d_tfp16); std(d_tfp16); corr(d_tfp16,d_tfp07)];
tfp_momentnames = {'mean'; 'stdev'; 'corr(.,d_tfp_07)'};
Table_dtfp_stats = table(d_tfp_07,d_tfp_11,d_tfp_13,d_tfp_14,d_tfp_15,d_tfp_16,'Rownames',tfp_momentnames)

d_u_07 = [mean(d_u07); std(d_u07); corr(d_u07,d_u07)];
d_u_11 = [mean(d_u11); std(d_u11); corr(d_u11,d_u07)];
d_u_13 = [mean(d_u13); std(d_u13); corr(d_u13,d_u07)];
d_u_14 = [mean(d_u14); std(d_u14); corr(d_u14,d_u07)];
d_u_15 = [mean(d_u15); std(d_u15); corr(d_u15,d_u07)];
d_u_16 = [mean(d_u16); std(d_u16); corr(d_u16,d_u07)];
u_momentnames = {'mean'; 'stdev'; 'corr(.,d_u_07)'};
Table_du_stats = table(d_u_07,d_u_11,d_u_13,d_u_14,d_u_15,d_u_16,'Rownames',u_momentnames)

dy = [mean(d_y); std(d_y); corr(d_y,d_y)];
dc = [mean(d_c); std(d_c); corr(d_c,d_y)];
di = [mean(d_I); std(d_I); corr(d_I,d_y)];
dH = [mean(d_h); std(d_h); corr(d_h,d_y)];
dh_n = [mean(d_h_n); std(d_h_n); corr(d_h_n,d_y)];
nipa_momentnames = {'mean'; 'stdev'; 'corr(.,d_y)'};
Table_nipa_stats = table(dy,dc,di,dH,dh_n,'Rownames',nipa_momentnames)

corr_nipa_07 = corr([d_y d_c d_I d_h d_h_n d_tfp07 d_u07 d_tfpu07]);
corr_nipa_11 = corr([d_y d_c d_I d_h d_h_n d_tfp11 d_u11 d_tfpu11]);
corr_nipa_13 = corr([d_y d_c d_I d_h d_h_n d_tfp13 d_u13 d_tfpu13]);
corr_nipa_14 = corr([d_y d_c d_I d_h d_h_n d_tfp14 d_u14 d_tfpu14]);
corr_nipa_15 = corr([d_y d_c d_I d_h d_h_n d_tfp15 d_u15 d_tfpu15]);
corr_nipa_16 = corr([d_y d_c d_I d_h d_h_n d_tfp16 d_u16 d_tfpu16]);

d_tfpu_07 = corr_nipa_07(:,end);
d_tfpu_11 = corr_nipa_11(:,end);
d_tfpu_13 = corr_nipa_13(:,end);
d_tfpu_14 = corr_nipa_14(:,end);
d_tfpu_15 = corr_nipa_15(:,end);
d_tfpu_16 = corr_nipa_16(:,end);
nipa_corrnames = {'corr(.,dy)'; 'corr(.,dc)'; 'corr(.,dI)'; 'corr(.,dh)'; 'corr(.,dh_n)'; 'corr(.,dtfp_)'; 'corr(.,du_)'; 'corr(.,dtfpu_)'};
Table_dtfpu_nipa = table(d_tfpu_07,d_tfpu_11,d_tfpu_13,d_tfpu_14,d_tfpu_15,d_tfpu_16,'Rownames',nipa_corrnames)

d_tfp_07 = corr_nipa_07(:,end-2);
d_tfp_11 = corr_nipa_11(:,end-2);
d_tfp_13 = corr_nipa_13(:,end-2);
d_tfp_14 = corr_nipa_14(:,end-2);
d_tfp_15 = corr_nipa_15(:,end-2);
d_tfp_16 = corr_nipa_16(:,end-2);
Table_dtfp_nipa = table(d_tfp_07,d_tfp_11,d_tfp_13,d_tfp_14,d_tfp_15,d_tfp_16,'Rownames',nipa_corrnames)

d_u_07 = corr_nipa_07(:,end-1);
d_u_11 = corr_nipa_11(:,end-1);
d_u_13 = corr_nipa_13(:,end-1);
d_u_14 = corr_nipa_14(:,end-1);
d_u_15 = corr_nipa_15(:,end-1);
d_u_16 = corr_nipa_16(:,end-1);
Table_du_nipa = table(d_u_07,d_u_11,d_u_13,d_u_14,d_u_15,d_u_16,'Rownames',nipa_corrnames)


%% create synthetic TFPs using different utilization series (for Table 1 in appendix)
disp('Unconditional moments for Table 1 of appendix')

sy_tfp07_07 = d_tfp07 - d_u07; 
sy_tfp11_07 = d_tfp11 - d_u07;
sy_tfp13_07 = d_tfp13 - d_u07;
sy_tfp14_07 = d_tfp14 - d_u07;
sy_tfp15_07 = d_tfp15 - d_u07;
sy_tfp16_07 = d_tfp16 - d_u07;

sy_tfp07_13 = d_tfp07 - d_u13; 
sy_tfp11_13 = d_tfp11 - d_u13;
sy_tfp13_13 = d_tfp13 - d_u13;
sy_tfp14_13 = d_tfp14 - d_u13;
sy_tfp15_13 = d_tfp15 - d_u13;
sy_tfp16_13 = d_tfp16 - d_u13;

sy_tfp07_14 = d_tfp07 - d_u14; 
sy_tfp11_14 = d_tfp11 - d_u14;
sy_tfp13_14 = d_tfp13 - d_u14;
sy_tfp14_14 = d_tfp14 - d_u14;
sy_tfp15_14 = d_tfp15 - d_u14;
sy_tfp16_14 = d_tfp16 - d_u14;

sy_tfp07_16 = d_tfp07 - d_u16; 
sy_tfp11_16 = d_tfp11 - d_u16;
sy_tfp13_16 = d_tfp13 - d_u16;
sy_tfp14_16 = d_tfp14 - d_u16;
sy_tfp15_16 = d_tfp15 - d_u16;
sy_tfp16_16 = d_tfp16 - d_u16;

sy_corr07 = corr([sy_tfp07_07 sy_tfp11_07 sy_tfp13_07 sy_tfp14_07 sy_tfp15_07 sy_tfp16_07]);
stfpu_07_07 = sy_corr07(:,1);
stfpu_11_07 = sy_corr07(:,2);
stfpu_13_07 = sy_corr07(:,3);
stfpu_14_07 = sy_corr07(:,4);
stfpu_15_07 = sy_corr07(:,5);
stfpu_16_07 = sy_corr07(:,6);
sy_corrnames = {'corr(.,stfpu07_07)'; 'corr(.,stfpu11_07)'; 'corr(.,stfpu13_07)'; 'corr(.,stfpu14_07)'; 'corr(.,stfpu15_07)'; 'corr(.,stfpu16_07)'};
sy_07_table = table(stfpu_07_07,stfpu_11_07,stfpu_13_07,stfpu_14_07,stfpu_15_07,stfpu_16_07,'Rownames',sy_corrnames)

sy_corr13 = corr([sy_tfp07_13 sy_tfp11_13 sy_tfp13_13 sy_tfp14_13 sy_tfp15_13 sy_tfp16_13]);
stfpu_07_13 = sy_corr13(:,1);
stfpu_11_13 = sy_corr13(:,2);
stfpu_13_13 = sy_corr13(:,3);
stfpu_14_13 = sy_corr13(:,4);
stfpu_15_13 = sy_corr13(:,5);
stfpu_16_13 = sy_corr13(:,6);
sy_corrnames = {'corr(.,stfpu07_13)'; 'corr(.,stfpu11_13)'; 'corr(.,stfpu13_13)'; 'corr(.,stfpu14_13)'; 'corr(.,stfpu15_13)'; 'corr(.,stfpu16_13)'};
sy_13_table = table(stfpu_07_13,stfpu_11_13,stfpu_13_13,stfpu_14_13,stfpu_15_13,stfpu_16_13,'Rownames',sy_corrnames)

sy_corr14 = corr([sy_tfp07_14 sy_tfp11_14 sy_tfp13_14 sy_tfp14_14 sy_tfp15_14 sy_tfp16_14]);
stfpu_07_14 = sy_corr14(:,1);
stfpu_11_14 = sy_corr14(:,2);
stfpu_13_14 = sy_corr14(:,3);
stfpu_14_14 = sy_corr14(:,4);
stfpu_15_14 = sy_corr14(:,5);
stfpu_16_14 = sy_corr14(:,6);
sy_corrnames = {'corr(.,stfpu07_14)'; 'corr(.,stfpu11_14)'; 'corr(.,stfpu13_14)'; 'corr(.,stfpu14_14)'; 'corr(.,stfpu15_14)'; 'corr(.,stfpu16_14)'};
sy_14_table = table(stfpu_07_14,stfpu_11_14,stfpu_13_14,stfpu_14_14,stfpu_15_14,stfpu_16_14,'Rownames',sy_corrnames)

sy_corr16 = corr([sy_tfp07_16 sy_tfp11_16 sy_tfp13_16 sy_tfp14_16 sy_tfp15_16 sy_tfp16_16]);
stfpu_07_16 = sy_corr16(:,1);
stfpu_11_16 = sy_corr16(:,2);
stfpu_13_16 = sy_corr16(:,3);
stfpu_14_16 = sy_corr16(:,4);
stfpu_15_16 = sy_corr16(:,5);
stfpu_16_16 = sy_corr16(:,6);
sy_corrnames = {'corr(.,stfpu07_16)'; 'corr(.,stfpu11_16)'; 'corr(.,stfpu13_16)'; 'corr(.,stfpu14_16)'; 'corr(.,stfpu15_16)'; 'corr(.,stfpu16_16)'};
sy_16_table = table(stfpu_07_16,stfpu_11_16,stfpu_13_16,stfpu_14_16,stfpu_15_16,stfpu_16_16,'Rownames',sy_corrnames)


%% Synethetic data for counterfactual TFPu and u series
tfp_u_synthetic = xlsread('KurmannSims_Part1.xlsx','tfp_util_synthetic','B5:J269');
u_synthetic = xlsread('KurmannSims_Part1.xlsx','util_synthetic','B5:J269');
[T,q] = size(tfp_u_data);
st = 1; %1 for 1947q3, 51 for 1960q1, 147 for 1984q1
et = 241; %T=265 for 2013q3, 130 for 1979q4, 241 for 2007q3
time = [1947.75:0.25:2007.75]; %adjust accordingly

d_u13_bw_bffk = u_synthetic(st:et,1);   %using bi-weight filtering and bffk industry weights and proportionality factors
d_u13_bw = u_synthetic(st:et,2);        %using bi-weight filtering and orginal bfk industry weights and proportionality factors
d_u13_bffk = u_synthetic(st:et,3);      %using original bp filtering and bffk industry weights and proportionality factors
d_u13_bpf_8_32 = u_synthetic(st:et,4);  %using bp(8,32) filtering instead of bp(8,160) which is what Fernald mistakenly uses for 2013.12 vintage
d_u13_nofilter = u_synthetic(st:et,5);  %using no filter
d_u13_bw2 = u_synthetic(st:et,6);       %using bi-weight filtering on hours level as opposed to hours differences (consistent with what Fernald does for bp filtered vintages)
u13_bw_bffk = cumsum(d_u13_bw_bffk/400);
u13_bw = cumsum(d_u13_bw/400);
u13_bffk = cumsum(d_u13_bffk/400);
u13_bpf_8_32 = cumsum(d_u13_bpf_8_32/400);
u13_nofilter = cumsum(d_u13_nofilter/400);
u13_bw2 = cumsum(d_u13_bw2/400);

d_tfpu13_bw_bffk = tfp_u_synthetic(st:et,1);
d_tfpu13_bw = tfp_u_synthetic(st:et,2);      
d_tfpu13_bffk = tfp_u_synthetic(st:et,3);
d_tfpu13_bpf_8_32 = tfp_u_synthetic(st:et,4);
d_tfpu13_nofilter = tfp_u_synthetic(st:et,5);
d_tfpu13_bw2 = tfp_u_synthetic(st:et,6);     
tfpu13_bw_bffk = cumsum(d_tfpu13_bw_bffk/400);
tfpu13_bw = cumsum(d_tfpu13_bw/400);
tfpu13_bffk = cumsum(d_tfpu13_bffk/400);
tfpu13_bpf_8_32 = cumsum(d_tfpu13_bpf_8_32/400);
tfpu13_nofilter = cumsum(d_tfpu13_nofilter/400);
tfpu13_bw2 = cumsum(d_tfpu13_bw2/400);

%% unconditional moments for Table 2 
disp('Unconditional moments for Table 2 of main text')

d_u_07 = [mean(d_u07); std(d_u07); corr(d_u07,d_u07)];
d_u_13 = [mean(d_u13); std(d_u13); corr(d_u13,d_u07)];
d_u_13_bpf_8_32 = [mean(d_u13_bpf_8_32); std(d_u13_bpf_8_32); corr(d_u13_bpf_8_32,d_u07)];
d_u_13_bw = [mean(d_u13_bw); std(d_u13_bw); corr(d_u13_bw,d_u07)];
d_u_13_bffk = [mean(d_u13_bffk); std(d_u13_bffk); corr(d_u13_bffk,d_u07)];
d_u_13_bw_bffk = [mean(d_u13_bw_bffk); std(d_u13_bw_bffk); corr(d_u13_bw_bffk,d_u07)];
d_u_14 = [mean(d_u14); std(d_u14); corr(d_u14,d_u07)];
d_u_13_nofilter = [mean(d_u13_nofilter); std(d_u13_nofilter); corr(d_u13_nofilter,d_u07)];
u_momentnames = {'mean'; 'stdev'; 'corr(.,d_u_07)'};
Table_du_stats = table(d_u_07,d_u_13,d_u_13_bffk,d_u_13_bw,d_u_13_bw_bffk,d_u_14,'Rownames',u_momentnames)
Table_du_stats_addtl = table(d_u_13,d_u_13_bpf_8_32,d_u_13_bw,d_u_13_nofilter,'Rownames',u_momentnames)

d_tfpu_07 = [mean(d_tfpu07); std(d_tfpu07); corr(d_tfpu07,d_tfpu07)];
d_tfpu_13 = [mean(d_tfpu13); std(d_tfpu13); corr(d_tfpu13,d_tfpu07)];
d_tfpu_13_bpf_8_32 = [mean(d_tfpu13_bpf_8_32); std(d_tfpu13_bpf_8_32); corr(d_tfpu13_bpf_8_32,d_tfpu07)];
d_tfpu_13_bw = [mean(d_tfpu13_bw); std(d_tfpu13_bw); corr(d_tfpu13_bw,d_tfpu07)];
d_tfpu_13_bffk = [mean(d_tfpu13_bffk); std(d_tfpu13_bffk); corr(d_tfpu13_bffk,d_tfpu07)];
d_tfpu_13_bw_bffk = [mean(d_tfpu13_bw_bffk); std(d_tfpu13_bw_bffk); corr(d_tfpu13_bw_bffk,d_tfpu07)];
d_tfpu_14 = [mean(d_tfpu14); std(d_tfpu14); corr(d_tfpu14,d_tfpu07)];
d_tfpu_13_nofilter = [mean(d_tfpu13_nofilter); std(d_tfpu13_nofilter); corr(d_tfpu13_nofilter,d_tfpu07)];
tfpu_momentnames = {'mean'; 'stdev'; 'corr(.,d_tfpu_07)'};
Table_dtfpu_stats = table(d_tfpu_07,d_tfpu_13,d_tfpu_13_bffk,d_tfpu_13_bw,d_tfpu_13_bw_bffk,d_tfpu_14,'Rownames',tfpu_momentnames)
Table_dtfpu_stats_addtl = table(d_tfpu_13,d_tfpu_13_bpf_8_32,d_tfpu_13_bw,d_tfpu_13_nofilter,'Rownames',tfpu_momentnames)
                
                
%% Figure 1 of Appendix plus additional figures to check results

% Figure 1 Appendix - comparison of 2013 vintage without filtering of hours and 2013 vintage bi-weight filtered
figure('Name','Figure 2 of Appendix')
subplot(211); fh = NBERbc(time, [u13_bw-mean(u13_bw)*ones(T,1) u13_nofilter-mean(u13_nofilter)*ones(T,1)], {'-','--'}, [1.5,1.5], {'k','r'}); 
                title('Utilization'); 
                legend([fh(1),fh(2)],'2013 vintage BW filtered','2013 vintage unfiltered')
subplot(212);   fh = NBERbc(time, [d_tfpu13_bw d_tfpu13_nofilter], {'-','--'}, [1.5,1.5], {'k','r'});
                title('Utilization-adjusted TFP growth');
                legend([fh(1),fh(2)],'2013 vintage BW filtered','2013 vintage unfiltered')
                
% Additional figures  

% Figure showing close match of 2013 vintage with bi-weight filtering,
% BFFK proportionality factors and industry weights relative to 2014 vintage
figure('Name','Figure 2_add1 for Appendix')
subplot(211);   fh = NBERbc(time, [u13-mean(u13)*ones(T,1) u13_bw_bffk-mean(u13_bw_bffk)*ones(T,1) u14-mean(u14)*ones(T,1)], {'-','-',':'}, [1.5,1.5,1.5], {'k','b','r'}); 
                title('Utilization'); 
                legend([fh(1),fh(2),fh(3)],'2013 vintage','2013 vintage BW BFFK','2014 vintage')
subplot(212);   fh = NBERbc(time, [d_tfpu13 d_tfpu13_bw_bffk d_tfpu14], {'-','-',':'}, [1.5,1.5,1.5], {'k','b','r'}); 
                title('Utilization-adjusted TFP growth');
                legend([fh(1),fh(2),fh(3)],'2013 vintage','2013 vintage BW BFFK','2014 vintage')

% Figure showing close match of 2013 vintage with 2013 vintage using BFFK proportionality factors and industry weights
figure('Name','Figure 2_add2 for Appendix')
subplot(211);   fh = NBERbc(time, [u13-mean(u13)*ones(T,1) u13_bffk-mean(u13_bffk)*ones(T,1)], {'-','--'}, [1.5,1.5], {'k','r'}); 
                title('Utilization'); 
                legend([fh(1),fh(2)],'2013 vintage','2013 vintage BFFK')
subplot(212);   fh = NBERbc(time, [d_tfpu13-mean(d_tfpu13)*ones(T,1) d_tfpu13_bffk-mean(d_tfpu13_bffk)*ones(T,1)], {'-','--'}, [1.5,1.5], {'k','r'}); 
                title('Utilization-adjusted TFP growth'); 
                legend([fh(1),fh(2)],'2013 vintage','2013 vintage BFFK')
  
                
%% Figure 2 for Appendix 

%detrending log average weekly hours in two different ways
bwtrend_h_n = bw_trend(h_n,48);      %trend from bi-weight filter that Fernald uses based on Stock-Watson's original code
cf_h_n = cf_bpass_default(h_n,8,32); %CF filtered series
cftrend_h_n = h_n - cf_h_n;          %trend implied by CF filter
bw_h_n = h_n - bwtrend_h_n;          %BW filtered series

figure('Name','Figure 3 for Appendix')
subplot(211);   fh = NBERbc(time, [h_n cftrend_h_n bwtrend_h_n], {'-','-','-'}, [1.5,1.5,1.5], {'r','k','b'}); 
                axis([time(1) time(end) 4.6 4.8])
                title('Log weekly hours');
                legend([fh(1),fh(2),fh(3)],'log weekly hours','CF 8-32 trend','BW 48 trend')
subplot(212);   fh = NBERbc(time, [cf_h_n bw_h_n], {'-','-'}, [1.5,1.5], {'k','b'}); 
                title('CF-filtered and BW-filtered log weekly hours')
                legend([fh(1),fh(2)],'CF 8-32 filtered','BW 48 filtered')
        
             
                

                











                
                
                
